The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file
Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here
#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))
ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))
processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd
|
| | 0%
|
|.... | 4%
|
|......... | 9%
|
|............. | 13%
|
|.................. | 17%
|
|...................... | 22%
|
|.......................... | 26%
|
|............................... | 30%
|
|................................... | 35%
|
|........................................ | 39%
|
|............................................ | 43%
|
|................................................ | 48%
|
|..................................................... | 52%
|
|......................................................... | 57%
|
|............................................................. | 61%
|
|.................................................................. | 65%
|
|...................................................................... | 70%
|
|........................................................................... | 74%
|
|............................................................................... | 78%
|
|................................................................................... | 83%
|
|........................................................................................ | 87%
|
|............................................................................................ | 91%
|
|................................................................................................. | 96%
|
|.....................................................................................................| 100%
output file: /tmp/RtmpswludR/file3e673faf9754
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Attaching package: ‘flextable’
The following object is masked from ‘package:purrr’:
compose
Attaching package: ‘ftExtra’
The following object is masked from ‘package:flextable’:
separate_header
Warning: Assuming taxa are rows
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3
library(cowplot)
library(flextable)
library(ftExtra)
This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.
Reshape back into an ASV matrix, but this time correcting for total abundance
raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")
raMat1 <- raMat %>% as.matrix()
countMat <- nonSpikes_Remake %>%
pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
column_to_rownames("ID") %>% as.matrix()
min(seqDep)
[1] 852
This value is lower than the lowist chimera checked value because the spikes have been discarded (while chimera checked read depth still has spikes)
sampleRichness <- rarefy(countMat, min(seqDep))
rarefy everything to the minimum depth (852)
countRare <- rrarefy(countMat, min(seqDep))
Gamma diversity
specpool(countRare)
Doesn’t finish
#specpool(countMat)
All richness estimates
richnessRare <- estimateR(countRare)
Shannon diversity
shan <- diversity(countRare)
shan
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2 3-1-S-180
4.447712 5.116637 4.678407 5.866088 5.097057 3.871649 5.479328 4.615634 4.929580 4.704532
3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500 3-2-B-53 3-2-S-0-2
5.260828 4.837697 4.502874 4.691918 4.620095 5.164912 5.416675 4.988461 4.938664 3.732937
3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2 3-3-B-1-2 3-3-B-180 3-3-B-20
4.849386 4.751104 4.843166 5.166633 4.808113 4.307320 4.350618 4.971685 3.437542 5.685777
3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500 3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180
5.260559 5.262516 5.343435 4.981628 4.942967 4.859932 4.357677 4.296520 5.024896 4.583716
4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2 4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2
4.394405 4.630198 4.426126 4.227551 4.968388 4.715856 5.195087 4.034619 4.583134 2.805419
4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53 5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53
4.563265 4.627756 4.805569 4.460398 4.284061 4.554498 4.210847 3.956005 4.066923 4.196692
5-5-B-0-2 5-5-B-180 5-5-B-5 5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2
4.573474 5.192896 5.419692 5.098760 4.829526 4.316809 5.003296 4.843684 4.216447 4.377762
C_5P1B_180 C_5P1B_1P2 C_5P1B_20 C_5P1B_500 C_5P1B_53
4.814288 4.856735 5.376345 4.676670 5.054771
Evenness
pielouJ <- shan/richnessRare["S.chao1",]
pielouJ
3-1-B-0-2 3-1-B-1-2 3-1-B-180 3-1-B-20 3-1-B-5 3-1-B-500 3-1-B-53 3-1-S-0-2 3-1-S-1-2
0.012242431 0.007628407 0.010461256 0.002668322 0.006980389 0.042083146 0.007982186 0.008122541 0.007524956
3-1-S-180 3-1-S-20 3-1-S-5 3-2-B-0-2 3-2-B-1-2 3-2-B-180 3-2-B-20 3-2-B-5 3-2-B-500
0.007157256 0.008075588 0.006058481 0.012802536 0.008688737 0.013388595 0.005490156 0.007066346 0.009079896
3-2-B-53 3-2-S-0-2 3-2-S-1-2 3-2-S-180 3-2-S-20 3-2-S-5 3-2-S-500 3-2-S-53 3-3-B-0-2
0.010109855 0.026272434 0.007768416 0.009249390 0.009327234 0.007932156 0.006879380 0.014943315 0.011706652
3-3-B-1-2 3-3-B-180 3-3-B-20 3-3-B-5 3-3-B-500 3-3-B-53 3-3-S-180 3-3-S-20 3-3-S-500
0.007745782 0.061937694 0.005071165 0.006867398 0.006900477 0.005323821 0.009684942 0.008590692 0.008665187
3-3-S-53 4-3-B-0-2 4-3-B-1-2 4-3-B-180 4-3-B-20 4-3-B-5 4-3-B-500 4-3-B-53 4-3-O-1-2
0.009491849 0.009200256 0.007283966 0.010619160 0.008898709 0.009710793 0.010432976 0.008489059 0.008707836
4-3-O-180 4-3-O-5 4-3-O-500 4-3-O-53 4-3-S-0-2 4-3-S-180 4-3-S-20 4-3-S-500 4-3-S-53
0.012332810 0.005953768 0.016115035 0.012607328 0.124685269 0.012810964 0.008107808 0.010446890 0.013387723
5-1-S-1-2 5-1-S-180 5-1-S-20 5-1-S-5 5-1-S-500 5-1-S-53 5-5-B-0-2 5-5-B-180 5-5-B-5
0.008619198 0.012719017 0.014697546 0.015157109 0.015985286 0.015201275 0.007584534 0.007587199 0.006467270
5-5-B-500 5-5-B-53 5-5-S-180 5-5-S-5 5-5-S-500 5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2
0.009734807 0.009680381 0.014795877 0.009898759 0.008377902 0.016317124 0.008321046 0.007915634 0.005932011
C_5P1B_20 C_5P1B_500 C_5P1B_53
0.003437984 0.008991717 0.006222376
diversityData <- sampleData %>%
left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
arrange(Size_Class)
Parameters for all plots
plotSpecs <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
Observed species counts, on rarefied data
plotObs <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs
Estemated Chao1 Richness
plotChao1 <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) +
scale_y_log10() +
ylab("Richness (Chao1)")
plotChao1
Shannon diversity
plotShan <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
ylab("Diversity (Shannon H)") +
lims(y = c(2.5, 6))
plotShan
Evenness
plotPielou <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou
All plots together
plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha
ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)
Rarefied observed species numbers
obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)
Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-221.223 -35.304 4.193 45.989 200.602
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 253678.231 131330.970 1.932 0.057517 .
log(Size_Class) 31.672 7.957 3.980 0.000168 ***
I(log(Size_Class)^2) -6.197 1.480 -4.186 8.24e-05 ***
lat -13212.822 6843.487 -1.931 0.057628 .
I(lat^2) 172.124 89.089 1.932 0.057461 .
depth 4.111 3.151 1.305 0.196386
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 74.8 on 69 degrees of freedom
Multiple R-squared: 0.2634, Adjusted R-squared: 0.21
F-statistic: 4.935 on 5 and 69 DF, p-value: 0.0006408
Rarified chao1 estimates
chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) +
I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-544.01 -152.34 -28.53 144.59 1401.31
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.037e+06 5.008e+05 2.071 0.04211 *
log(Size_Class) 8.827e+01 3.034e+01 2.909 0.00487 **
I(log(Size_Class)^2) -1.931e+01 5.645e+00 -3.421 0.00105 **
lat -5.406e+04 2.610e+04 -2.072 0.04204 *
I(lat^2) 7.044e+02 3.397e+02 2.074 0.04185 *
depth 2.076e+01 1.202e+01 1.727 0.08857 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 285.2 on 69 degrees of freedom
Multiple R-squared: 0.2082, Adjusted R-squared: 0.1509
F-statistic: 3.63 on 5 and 69 DF, p-value: 0.005664
As above but without latitude and depth
chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)
Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2),
data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-469.09 -166.48 -32.74 134.06 1512.72
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 597.962 51.919 11.517 < 2e-16 ***
log(Size_Class) 88.481 30.833 2.870 0.005391 **
I(log(Size_Class)^2) -19.759 5.736 -3.445 0.000956 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 290.2 on 72 degrees of freedom
Multiple R-squared: 0.1447, Adjusted R-squared: 0.121
F-statistic: 6.092 on 2 and 72 DF, p-value: 0.003596
shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(shanMod)
Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-1.33049 -0.22559 0.06304 0.28864 0.70832
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.424e+03 7.757e+02 1.836 0.070648 .
log(Size_Class) 1.914e-01 4.700e-02 4.072 0.000122 ***
I(log(Size_Class)^2) -3.501e-02 8.744e-03 -4.004 0.000155 ***
lat -7.395e+01 4.042e+01 -1.830 0.071627 .
I(lat^2) 9.627e-01 5.262e-01 1.830 0.071630 .
depth 1.458e-02 1.861e-02 0.783 0.436226
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4418 on 69 degrees of freedom
Multiple R-squared: 0.2814, Adjusted R-squared: 0.2294
F-statistic: 5.405 on 5 and 69 DF, p-value: 0.0002982
pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)
Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityData)
Residuals:
Min 1Q Median 3Q Max
-0.016051 -0.004926 -0.002385 0.001354 0.099434
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.747e+01 2.640e+01 -1.041 0.3016
log(Size_Class) -4.045e-03 1.599e-03 -2.529 0.0137 *
I(log(Size_Class)^2) 7.144e-04 2.976e-04 2.401 0.0191 *
lat 1.431e+00 1.375e+00 1.040 0.3019
I(lat^2) -1.861e-02 1.791e-02 -1.039 0.3023
depth -4.148e-04 6.334e-04 -0.655 0.5148
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01504 on 69 degrees of freedom
Multiple R-squared: 0.1054, Adjusted R-squared: 0.0406
F-statistic: 1.626 on 5 and 69 DF, p-value: 0.1646
uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2
predict(obsMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
230.1017 230.1017 201.2125 201.2125 209.8268 162.3905 162.3905 187.5034 200.5021 302.6958 302.6958 273.8066
13 14 15 16 17 18 19 20 21 22 23 24
273.8066 282.4209 234.9846 234.9846 260.0975 260.0975 332.0508 332.0508 303.1616 303.1616 311.7759 264.3396
25 26 27 28 29 30 31 32 33 34 35 36
264.3396 289.4525 302.4512 302.4512 336.3980 336.3980 307.5088 307.5088 316.1232 316.1232 268.6868 268.6868
37 38 39 40 41 42 43 44 45 46 47 48
293.7998 293.7998 325.1970 296.3078 296.3078 304.9221 304.9221 257.4858 257.4858 257.4858 282.5987 282.5987
49 50 51 52 53 54 55 56 57 58 59 60
295.5974 295.5974 294.4974 294.4974 265.6083 265.6083 274.2226 274.2226 226.7862 226.7862 226.7862 251.8992
61 62 63 64 65 66 67 68 69 70 71 72
251.8992 264.8979 264.8979 254.6368 225.7477 225.7477 234.3620 234.3620 186.9256 186.9256 186.9256 212.0386
73 74 75
212.0386 225.0373 225.0373
$se.fit
[1] 26.16949 26.16949 25.45599 25.45599 25.32015 26.93768 26.93768 28.45024 33.00461 18.35754 18.35754
[12] 17.75988 17.75988 16.67860 19.43184 19.43184 20.85030 20.85030 18.66951 18.66951 18.22171 18.22171
[23] 16.71287 19.47713 19.47713 20.56141 27.32163 27.32163 18.66740 18.66740 18.18962 18.18962 16.44692
[34] 16.44692 19.01561 19.01561 19.95080 19.95080 17.99853 17.37889 17.37889 15.51257 15.51257 17.88327
[45] 17.88327 17.88327 18.84420 18.84420 25.74850 25.74850 18.72596 18.72596 17.86591 17.86591 16.17413
[56] 16.17413 17.86737 17.86737 17.86737 18.90703 18.90703 25.32224 25.32224 23.71732 22.79314 22.79314
[67] 21.66683 21.66683 22.43854 22.43854 22.43854 23.40764 23.40764 28.33970 28.33970
$df
[1] 69
$residual.scale
[1] 74.80449
diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit
plotSpecs2 <- list(
facet_wrap(~Depth, ncol = 1) ,
theme_bw(base_size = 16) ,
#geom_point(size = 4) ,
geom_path(aes(color = as.factor(Station))) ,
scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
#scale_y_log10nice() ,
scale_shape_manual(values = rep(21:25, 2)) ,
scale_fill_viridis_d(option = "plasma") ,
scale_color_viridis_d(option = "plasma") ,
labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
theme(legend.position = "none",
plot.margin = unit(c(0, 0, 0, 0), "cm"),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = .5),
axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)
plotObs_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted ASVs")
plotObs_pred
predict(chao1Mod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
513.9062 513.9062 390.1285 390.1285 469.8081 281.2073 281.2073 408.9487 387.9514 721.4337 721.4337 597.6560
13 14 15 16 17 18 19 20 21 22 23 24
597.6560 677.3356 488.7348 488.7348 616.4761 616.4761 798.0222 798.0222 674.2445 674.2445 753.9241 565.3233
25 26 27 28 29 30 31 32 33 34 35 36
565.3233 693.0646 672.0674 672.0674 797.1065 797.1065 673.3289 673.3289 753.0084 753.0084 564.4077 564.4077
37 38 39 40 41 42 43 44 45 46 47 48
692.1490 692.1490 752.0344 628.2567 628.2567 707.9363 707.9363 519.3355 519.3355 519.3355 647.0769 647.0769
49 50 51 52 53 54 55 56 57 58 59 60
626.0796 626.0796 643.6118 643.6118 519.8342 519.8342 599.5137 599.5137 410.9130 410.9130 410.9130 538.6543
61 62 63 64 65 66 67 68 69 70 71 72
538.6543 517.6571 517.6571 508.7384 384.9608 384.9608 464.6403 464.6403 276.0396 276.0396 276.0396 403.7809
73 74 75
403.7809 382.7836 382.7836
$se.fit
[1] 99.78989 99.78989 97.06920 97.06920 96.55117 102.71918 102.71918 108.48688 125.85368 70.00126
[11] 70.00126 67.72225 67.72225 63.59909 74.09781 74.09781 79.50669 79.50669 71.19088 71.19088
[21] 69.48331 69.48331 63.72976 74.27049 74.27049 78.40508 104.18328 104.18328 71.18284 71.18284
[31] 69.36094 69.36094 62.71564 62.71564 72.51061 72.51061 76.07670 76.07670 68.63229 66.26945
[41] 66.26945 59.15278 59.15278 68.19276 68.19276 68.19276 71.85699 71.85699 98.18459 98.18459
[51] 71.40613 71.40613 68.12657 68.12657 61.67546 61.67546 68.13214 68.13214 68.13214 72.09658
[61] 72.09658 96.55916 96.55916 90.43926 86.91517 86.91517 82.62028 82.62028 85.56300 85.56300
[71] 85.56300 89.25838 89.25838 108.06537 108.06537
$df
[1] 69
$residual.scale
[1] 285.2456
diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit
plotChao1_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred
predict(shanMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9 10 11 12
4.510911 4.510911 4.363021 4.363021 4.319694 4.051892 4.051892 4.143951 4.388759 4.943313 4.943313 4.795423
13 14 15 16 17 18 19 20 21 22 23 24
4.795423 4.752096 4.484294 4.484294 4.576353 4.576353 5.126894 5.126894 4.979003 4.979003 4.935677 4.667874
25 26 27 28 29 30 31 32 33 34 35 36
4.667874 4.759933 5.004741 5.004741 5.168679 5.168679 5.020788 5.020788 4.977462 4.977462 4.709659 4.709659
37 38 39 40 41 42 43 44 45 46 47 48
4.801718 4.801718 5.117504 4.969614 4.969614 4.926287 4.926287 4.658485 4.658485 4.658485 4.750544 4.750544
49 50 51 52 53 54 55 56 57 58 59 60
4.995352 4.995352 4.959250 4.959250 4.811359 4.811359 4.768033 4.768033 4.500230 4.500230 4.500230 4.592289
61 62 63 64 65 66 67 68 69 70 71 72
4.592289 4.837097 4.837097 4.746740 4.598849 4.598849 4.555523 4.555523 4.287720 4.287720 4.287720 4.379779
73 74 75
4.379779 4.624587 4.624587
$se.fit
[1] 0.15456512 0.15456512 0.15035102 0.15035102 0.14954865 0.15910231 0.15910231 0.16803594 0.19493546
[10] 0.10842534 0.10842534 0.10489537 0.10489537 0.09850898 0.11477052 0.11477052 0.12314836 0.12314836
[19] 0.11026796 0.11026796 0.10762308 0.10762308 0.09871139 0.11503798 0.11503798 0.12144207 0.16137006
[28] 0.16137006 0.11025550 0.11025550 0.10743355 0.10743355 0.09714060 0.09714060 0.11231209 0.11231209
[37] 0.11783563 0.11783563 0.10630494 0.10264512 0.10264512 0.09162207 0.09162207 0.10562415 0.10562415
[46] 0.10562415 0.11129969 0.11129969 0.15207866 0.15207866 0.11060136 0.11060136 0.10552162 0.10552162
[55] 0.09552947 0.09552947 0.10553025 0.10553025 0.10553025 0.11167080 0.11167080 0.14956102 0.14956102
[64] 0.14008188 0.13462339 0.13462339 0.12797101 0.12797101 0.13252901 0.13252901 0.13252901 0.13825280
[73] 0.13825280 0.16738306 0.16738306
$df
[1] 69
$residual.scale
[1] 0.4418186
diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit
plotShannonH_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5) +
plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred
predict(pielouMod, se.fit = TRUE)
$fit
1 2 3 4 5 6 7 8 9
0.019088298 0.019088298 0.021985268 0.021985268 0.021029079 0.025251377 0.025251377 0.022216580 0.018901598
10 11 12 13 14 15 16 17 18
0.010013570 0.010013570 0.012910541 0.012910541 0.011954352 0.016176649 0.016176649 0.013141852 0.013141852
19 20 21 22 23 24 25 26 27
0.006067283 0.006067283 0.008964253 0.008964253 0.008008064 0.012230362 0.012230362 0.009195565 0.005880583
28 29 30 31 32 33 34 35 36
0.005880583 0.005020093 0.005020093 0.007917064 0.007917064 0.006960875 0.006960875 0.011183172 0.011183172
37 38 39 40 41 42 43 44 45
0.008148375 0.008148375 0.005927552 0.008824522 0.008824522 0.007868333 0.007868333 0.012090631 0.012090631
46 47 48 49 50 51 52 53 54
0.012090631 0.009055834 0.009055834 0.005740852 0.005740852 0.008985179 0.008985179 0.011882150 0.011882150
55 56 57 58 59 60 61 62 63
0.010925961 0.010925961 0.015148258 0.015148258 0.015148258 0.012113461 0.012113461 0.008798480 0.008798480
64 65 66 67 68 69 70 71 72
0.013178110 0.016075081 0.016075081 0.015118892 0.015118892 0.019341189 0.019341189 0.019341189 0.016306392
73 74 75
0.016306392 0.012991411 0.012991411
$se.fit
[1] 0.005259890 0.005259890 0.005116483 0.005116483 0.005089178 0.005414292 0.005414292 0.005718306
[9] 0.006633703 0.003689742 0.003689742 0.003569616 0.003569616 0.003352286 0.003905670 0.003905670
[17] 0.004190770 0.004190770 0.003752447 0.003752447 0.003662441 0.003662441 0.003359174 0.003914772
[25] 0.003914772 0.004132704 0.005491464 0.005491464 0.003752023 0.003752023 0.003655991 0.003655991
[33] 0.003305719 0.003305719 0.003822009 0.003822009 0.004009976 0.004009976 0.003617584 0.003493039
[41] 0.003493039 0.003117922 0.003117922 0.003594417 0.003594417 0.003594417 0.003787557 0.003787557
[49] 0.005175275 0.005175275 0.003763792 0.003763792 0.003590927 0.003590927 0.003250892 0.003250892
[57] 0.003591221 0.003591221 0.003591221 0.003800186 0.003800186 0.005089600 0.005089600 0.004767022
[65] 0.004581268 0.004581268 0.004354886 0.004354886 0.004509996 0.004509996 0.004509996 0.004704778
[73] 0.004704778 0.005696088 0.005696088
$df
[1] 69
$residual.scale
[1] 0.0150352
diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit
plot_pielouJ_pred <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred
plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
plotPredictions
ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)
plot_grid(plotObs, plotChao1, plotShan, plotPielou,
plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
nrow = 2, labels = LETTERS)
Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).
alphaSummary <- tibble(
metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(obsMod, chao1Mod, shanMod, pielouMod)
)
alphaSummary <- alphaSummary %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
alphaSummary <- alphaSummary %>%
select(-model) %>%
unnest(df)
alphaSummary <- alphaSummary %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
bold(i = ~ p< 0.05, j = "p") %>%
colformat_md() %>%
set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
Metric | Term | Estimate | Standard | T Value | p |
Observed ASVs | Intercept | 2.5 x 105 | 1.3 x 105 | 1.93 | 0.058 |
log(Size Class) | 3.2 x 101 | 8.0 x 100 | 3.98 | < 0.001 | |
log(Size Class)2 | -6.2 x 100 | 1.5 x 100 | -4.19 | < 0.001 | |
Latitude | -1.3 x 104 | 6.8 x 103 | -1.93 | 0.058 | |
Latitude2 | 1.7 x 102 | 8.9 x 101 | 1.93 | 0.057 | |
Depth | 4.1 x 100 | 3.2 x 100 | 1.30 | 0.196 | |
Richness (Chao1) | Intercept | 1.0 x 106 | 5.0 x 105 | 2.07 | 0.042 |
log(Size Class) | 8.8 x 101 | 3.0 x 101 | 2.91 | 0.005 | |
log(Size Class)2 | -1.9 x 101 | 5.6 x 100 | -3.42 | 0.001 | |
Latitude | -5.4 x 104 | 2.6 x 104 | -2.07 | 0.042 | |
Latitude2 | 7.0 x 102 | 3.4 x 102 | 2.07 | 0.042 | |
Depth | 2.1 x 101 | 1.2 x 101 | 1.73 | 0.089 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.8 x 102 | 1.84 | 0.071 |
log(Size Class) | 1.9 x 10-1 | 4.7 x 10-2 | 4.07 | < 0.001 | |
log(Size Class)2 | -3.5 x 10-2 | 8.7 x 10-3 | -4.00 | < 0.001 | |
Latitude | -7.4 x 101 | 4.0 x 101 | -1.83 | 0.072 | |
Latitude2 | 9.6 x 10-1 | 5.3 x 10-1 | 1.83 | 0.072 | |
Depth | 1.5 x 10-2 | 1.9 x 10-2 | 0.78 | 0.436 | |
Evenness (Pielou J) | Intercept | -2.7 x 101 | 2.6 x 101 | -1.04 | 0.302 |
log(Size Class) | -4.0 x 10-3 | 1.6 x 10-3 | -2.53 | 0.014 | |
log(Size Class)2 | 7.1 x 10-4 | 3.0 x 10-4 | 2.40 | 0.019 | |
Latitude | 1.4 x 100 | 1.4 x 100 | 1.04 | 0.302 | |
Latitude2 | -1.9 x 10-2 | 1.8 x 10-2 | -1.04 | 0.302 | |
Depth | -4.1 x 10-4 | 6.3 x 10-4 | -0.65 | 0.515 |
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))
Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.
diversityDataWB <- full_join(diversityData,
richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
mutate(pielouJ2 = shannonH/break_estimate) %>%
identity()
diversityDataWB
pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)
Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-0.013975 -0.005102 -0.002495 0.000889 0.105935
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.971e+01 2.752e+01 -0.716 0.4762
log(Size_Class) -3.292e-03 1.668e-03 -1.974 0.0524 .
I(log(Size_Class)^2) 5.752e-04 3.103e-04 1.854 0.0680 .
lat 1.025e+00 1.434e+00 0.715 0.4771
I(lat^2) -1.332e-02 1.867e-02 -0.713 0.4780
depth -2.399e-04 6.605e-04 -0.363 0.7176
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.01568 on 69 degrees of freedom
Multiple R-squared: 0.06743, Adjusted R-squared: -0.000147
F-statistic: 0.9978 on 5 and 69 DF, p-value: 0.4257
Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.
plotBreak <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Richness (Breakaway)")
plotBreak
plotPielou2 <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
#scale_y_log10()+
ylab("Evenness (PielouJ)")
plotPielou2
predict(pielouMod2, se.fit = TRUE)
$fit
1 2 3 4 5 6 7
0.0116808693 0.0116808693 0.0135879652 0.0135879652 0.0133789427 0.0160393203 0.0160393203
8 9 10 11 12 13 14
0.0139920245 0.0099168119 0.0043118512 0.0043118512 0.0062189471 0.0062189471 0.0060099245
15 16 17 18 19 20 21
0.0086703022 0.0086703022 0.0066230064 0.0066230064 0.0010848943 0.0010848943 0.0029919902
22 23 24 25 26 27 28
0.0029919902 0.0027829677 0.0054433453 0.0054433453 0.0033960495 -0.0006791631 -0.0006791631
29 30 31 32 33 34 35
0.0001937160 0.0001937160 0.0021008119 0.0021008119 0.0018917894 0.0018917894 0.0045521670
36 37 38 39 40 41 42
0.0045521670 0.0025048712 0.0025048712 0.0008906780 0.0027977739 0.0027977739 0.0025887513
43 44 45 46 47 48 49
0.0025887513 0.0052491290 0.0052491290 0.0052491290 0.0032018332 0.0032018332 -0.0008733794
50 51 52 53 54 55 56
-0.0008733794 0.0033103787 0.0033103787 0.0052174746 0.0052174746 0.0050084520 0.0050084520
57 58 59 60 61 62 63
0.0076688297 0.0076688297 0.0076688297 0.0056215339 0.0056215339 0.0015463213 0.0015463213
64 65 66 67 68 69 70
0.0066511872 0.0085582831 0.0085582831 0.0083492605 0.0083492605 0.0110096382 0.0110096382
71 72 73 74 75
0.0110096382 0.0089623424 0.0089623424 0.0048871298 0.0048871298
$se.fit
[1] 0.005484417 0.005484417 0.005334888 0.005334888 0.005306418 0.005645409 0.005645409 0.005962400
[9] 0.006916873 0.003847244 0.003847244 0.003721991 0.003721991 0.003495383 0.004072389 0.004072389
[17] 0.004369659 0.004369659 0.003912625 0.003912625 0.003818778 0.003818778 0.003502565 0.004081880
[25] 0.004081880 0.004309115 0.005725876 0.005725876 0.003912183 0.003912183 0.003812053 0.003812053
[33] 0.003446829 0.003446829 0.003985157 0.003985157 0.004181148 0.004181148 0.003772006 0.003642145
[41] 0.003642145 0.003251015 0.003251015 0.003747850 0.003747850 0.003747850 0.003949234 0.003949234
[49] 0.005396190 0.005396190 0.003924456 0.003924456 0.003744212 0.003744212 0.003389661 0.003389661
[57] 0.003744518 0.003744518 0.003744518 0.003962402 0.003962402 0.005306857 0.005306857 0.004970510
[65] 0.004776827 0.004776827 0.004540781 0.004540781 0.004702512 0.004702512 0.004702512 0.004905608
[73] 0.004905608 0.005939234 0.005939234
$df
[1] 69
$residual.scale
[1] 0.015677
diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit
plot_pielouJ2_pred <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred
plotBreakaway <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
plotSpecs +
geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) +
scale_y_log10() +
ylab("Richness (Breakaway)")
plotBreakaway
#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame
Why are these not smooth curves?!! What if I redo the model, this time with the same data frame
breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat + I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()
Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) +
lat + I(lat^2) + depth, data = diversityDataWB)
Residuals:
Min 1Q Median 3Q Max
-2974.5 -1191.2 -151.6 599.9 6768.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7124615.61 3339862.88 2.133 0.0365 *
log(Size_Class) 244.45 202.35 1.208 0.2312
I(log(Size_Class)^2) -75.16 37.65 -1.996 0.0498 *
lat -370568.38 174035.93 -2.129 0.0368 *
I(lat^2) 4817.28 2265.61 2.126 0.0371 *
depth 151.10 80.15 1.885 0.0636 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared: 0.1414, Adjusted R-squared: 0.0792
F-statistic: 2.273 on 5 and 69 DF, p-value: 0.0567
Note the non statistical significance overall
#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit
plot_break_pred <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
# filter(Station == 4.3) %>%
ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both")) +
plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred
plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB
ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlphaWB, width = 8, height = 4)
Summary table I want both breakaway metrics here
bettaTable <- myBet$table %>%
as.data.frame() %>%
rename(estimate = Estimates,
`std.error` = `Standard Errors`,
`p.value`=`p-values`
) %>%
mutate(`statistic` = NA) %>%
rownames_to_column(var = "term") %>%
select(term, estimate, std.error, statistic, p.value) %>%
as_tibble()
bettaTable
alphaSummary2 <- tibble(
metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
model = list(breakLm, shanMod, pielouMod2)
)
alphaSummary2 <- alphaSummary2 %>%
mutate(df = map(model, ~broom::tidy(summary(.))))
## Add in willis variables
breakawaySummary <- tibble(
metric = "Richness (Breakaway -- Betta)",
model = NULL,
df = list(bettaTable)
)
alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)
alphaSummary2 <- alphaSummary2 %>%
select(-model) %>%
unnest(df)
alphaSummary2 <- alphaSummary2 %>%
rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
Term = str_replace(Term, "\\^2", "\\^2\\^"),
Term = str_replace(Term, "depth", "Depth"),
Term = str_replace(Term, "lat", "Latitude"),
Term = str_replace(Term, "_", " ")# BOOKMARK!!
) %>%
mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
reformat_sci()
) %>%
mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
rename(`Standard\nError` = `Standard Error`) %>%
identity()
alphaSummary2
alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
align(j = -c(1:2), align = "right")
alphaTable2
Metric | Term | Estimate | Standard | T Value | p |
Richness (Breakaway – Betta) | Intercept | 7.1 x 106 | 2.4 x 102 | NA | < 0.001 |
log(Size Class) | 1.2 x 102 | 6.1 x 101 | NA | 0.058 | |
log(Size Class)2 | -5.0 x 101 | 1.2 x 101 | NA | < 0.001 | |
Latitude | -3.7 x 105 | 6.1 x 100 | NA | < 0.001 | |
Latitude2 | 4.8 x 103 | 1.6 x 10-1 | NA | < 0.001 | |
Depth | 1.5 x 102 | 1.0 x 101 | NA | < 0.001 | |
Richness (Breakaway – LM) | Intercept | 7.1 x 106 | 3.3 x 106 | 2.13 | 0.036 |
log(Size Class) | 2.4 x 102 | 2.0 x 102 | 1.21 | 0.231 | |
log(Size Class)2 | -7.5 x 101 | 3.8 x 101 | -2.00 | 0.050 | |
Latitude | -3.7 x 105 | 1.7 x 105 | -2.13 | 0.037 | |
Latitude2 | 4.8 x 103 | 2.3 x 103 | 2.13 | 0.037 | |
Depth | 1.5 x 102 | 8.0 x 101 | 1.89 | 0.064 | |
Diversity (Shannon H) | Intercept | 1.4 x 103 | 7.8 x 102 | 1.84 | 0.071 |
log(Size Class) | 1.9 x 10-1 | 4.7 x 10-2 | 4.07 | < 0.001 | |
log(Size Class)2 | -3.5 x 10-2 | 8.7 x 10-3 | -4.00 | < 0.001 | |
Latitude | -7.4 x 101 | 4.0 x 101 | -1.83 | 0.072 | |
Latitude2 | 9.6 x 10-1 | 5.3 x 10-1 | 1.83 | 0.072 | |
Depth | 1.5 x 10-2 | 1.9 x 10-2 | 0.78 | 0.436 | |
Evenness (Pielou J) | Intercept | -2.0 x 101 | 2.8 x 101 | -0.72 | 0.476 |
log(Size Class) | -3.3 x 10-3 | 1.7 x 10-3 | -1.97 | 0.052 | |
log(Size Class)2 | 5.8 x 10-4 | 3.1 x 10-4 | 1.85 | 0.068 | |
Latitude | 1.0 x 100 | 1.4 x 100 | 0.71 | 0.477 | |
Latitude2 | -1.3 x 10-2 | 1.9 x 10-2 | -0.71 | 0.478 | |
Depth | -2.4 x 10-4 | 6.6 x 10-4 | -0.36 | 0.718 |
alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))
myBet$table
plotAlphaWB_pred <- plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)
plotAlphaWB_pred
ggsave(here::here("Figures", "BreakawayAlphaPredictions.png"), plot = plotAlphaWB_pred, width = 8, height = 4)